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We use quantum Monte Carlo simulations to study the effect of disorder, in the form of a disordered 
chemical potential, on the phase diagram of the hard core bosonic Hubbard model in two dimensions. 
We find numerical evidence that in two dimensions, no matter how weak the disorder, it will always 
destroy the long range density wave order (checkerboard solid) present at half filling and strong 
nearest neighbor repulsion and replace it with a bose glass phase. We study the properties of this 
glassy phase including the superfluid density, energy gaps and the full Green's function. We also 
study the possibility of other localized phases at weak nearest neighbor repulsion, i.e. Anderson 
localization. We find that such a phase does not truly exist: The disorder must exceed a threshold 
before the bosons (at weak nn repulsion) are localized. The phase diagram for hard core bosons 
with disorder cannot be obtained easily from the soft core phase diagram discussed in the literature. 

PACS numbers: 74.76.-w, 74.40.+k, 73.43.Nq 



I. INTRODUCTION 

The two dimensional bosonic Hubbard model has been 
the subject of intense interest these past years because it 
is thought to capture many of the important qualitative 
features of two dimensional superconductors and super- 
fluids at very low temperature. For example, Helium 
atoms adsorbed on a surfacetl can clearly be described 
by bosons moving in a two dimensional environment. It 
is then natural to examine the role of disorder in local- 
izing the bosons and producing exotic phases such as a 
bose glass or a normal fluid at zero temperature. In the 
case of soft core bosons with contact repulsion, tha.bose 
glass phase was predicted and stutUpd theoreticallycl and 
subsequently verified numericallju'cl. 

Another reason for the increased interest in dis- 
ordered bosonic systems is a set of fascinating ex- 
periments on the superconducting-insulating transition 
suggesting the possibility p©£--a--janiversal conductance 
right at the transitionHQCl'B'BlljO. Several ideas, based 
on disordered bosonic Hubbard models, have been 
suggestedH ±Q-Mqplain these results. Extensive numerical 
simulationadOEj appear to support these ideas qualita- 
tively, although the numerical values of the conductance 
are not in agreement. 

The question of existence of a normal conducting state 
at zero temperature has regained momentum with re- 
cent experimental discoveriesll3. Attempts to explain this 
phase pmceed via models of disordered bosons, see for 
exampleESIlZl and references therein. 

Yet another reason to study bosons in external po- 
tentials (random or otherwise) are the recent fascinating 



experimentS|J3n atomic Bose-Einstein condensates on op- 
tical latticests. In many cases, such as this one, the rele- 
vant bosonic Hubbard model is the soft core one in others 
it is the hard core that is of interest. It is therefore in- 
teresting and important to expose and understand some 
of the important differences between these two cases. 

The paper is organized as follows. In Sec. II we will 
first present the hard core boson Hubbard model, our 
simulation algorithm and measurements. Then, in Sec. 



[II we will first review the phase diagram of the clean 
model before presenting our results on the disordered 
model at strong and weak near neighbor repulsions. Con- 
clusions and comments are in section IV. 



II. THE BOSON HUBBARD MODEL 

The hard core boson Hubbard Hamiltonian is given by 
H = — < ^^(fl; aj + aj fli) — /iifti 

(iJ) i 

+ ^i^ninj (1) 

(i,j> 

ai (flj ) are destruction (creation) operators of hard-core 
bosons on site i of a two dimensional square lattice, and 
rii is the boson number at site i while /Zi is the site depen- 
dent chemical potential. This is, therefore, a site depen- 
dent energy which models the disorder in the system. In 
the absence of disorder, becomes the normal chemical 
potential, /x. There are other ways to model disorder. 
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For example, one can have bond dependent hopping pa- 
rameter or near-neighbor interaction (Vi ij). It is 
thought, though not fuUy demonstrated, that these pos- 
sibihties fall in the same universality class. The hopping 
parameter is chosen to be t = 1 to fix the energy scale. 
V is the near neighbor interaction. 

To characterize the different phases, we need to mea- 
sure several physical quantities. A superfluid phase is 
characterized by the absence of long range density or- 
der and a non-vanishing superfluid (SF) density. The SF 
density, Ps is given by 

Ps = {W^)/2tp, (2) 

where W is the winding number of the phase of the bosM, 
wave function in one of the two spatial dimensionsB'ta 
and f3 — 1/kT. Long range density order (such as in 
the checkerboard solid) is characterized by the density- 
density correlation function, c(l), and the structure fac- 
tor, S'(q), its Fourier transform. They are given by 

c(l) = 

5(q) = ^^e^'i'ca), (3) 
1 

where rij is the occupancy at site j. In the presence of 
long range order, S{q) will diverge with the system size 
for a given ordering momentum, q*, which characterizes 
the ordered phase. For example, for checkerboard order, 
q* = (7r,7r). 

Two other very useful quantities are the equal time 
Green's function 

G(|j-i|) = (ajat), (4) 
and the Green's function in imaginary time, 

G(r) = (ai,.4o)- (5) 

In the superfluid phase, G(|j — i|) saturates at a nonzero 
value for large separations, while G{t) tends to zero expo- 
nentially thus yielding the quasiparticle excitation energy 
spectrum. In the simulations, G{\i — i|) was measured 
along the lattice axes. 

In the presence of disorder, we need to average over re- 
alizations of disorder in addition to the usual statistical 
average for a given realization. The number of realiza- 
tions we used depended on the size of the system but is 
typically a few hundred. We do our simulations using the 
stochastic series expansion (SSE) algorithm with worm 
updatescj. This algorithm is numerically exact without 
any discretization error. In addition it uses non-local 
updates, hence even large systems can be sampled effi- 
ciently. 

III. THE PHASE DIAGRAM 

The phase diagram of the bosonic Hubbard model with 
finite contact repulsion (i.e. soft core) and no nearest 
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FIG. 1: The phase diagram of the hard core bosonic Hubbard 
model in the absence of disorder (8 x 8, /3 = 14). Dashed lines 
indicate first order transitions, continuous lines second order 
transitions. 

neighbor (nn) repulsion has been studied extensively in 
one and two dimensions both with and without disorder. 
In the absence of disorder, for incommensurate particle 
fillings, one always has a superfluid phase. This phase 
disappears at commensurate fiUinss [p — 1,2,3...) when 
the onsite repulsion is large enoughcl. The resulting phase 
is an incompressilpJ|e-|Mott insulating phase which takes 
the form of lobesoaO in the (t/Vb, ^/Vb) plane, where Vb 
is the onsite repulsion. This phase is gapped, there is a 
substantial energy cost (increase in the chemical poten- 
tial) for adding a particle onto a commensurate phase. 
It was arguedH that in the presence of any amount of 
disorder, a new compressible insulating (i.e. localized) 
phase, the bose glass, is produced at incommensurate 
fillings and that for strong enough disorder the gapped 
phase disappears 9}rti*^feini'^^i^ subsequently con- 
firmed numericallyu'QEJ'ESO. 

The picture changes for hard core bosons with near 
neighbor repulsion, V. The bosons still form a super- 
fluid for incommensurate particle flUing and V not too 
strong. Also, at full filling, the bosons are always frozen 
into a Mott insulator since hopping to a neighbor would 
produce double occupancy which is strictly forbidden. 
At half filling, increasing V eventually freezes the bosons 
into an incompressible gapped checkerboard solid: Alter- 
nate sites are occupied since the presence of a neighbor 
costs too much energy and there is a big energy cost (gap) 
to add a particle. The phase diagram is shown in Fig. l|. 



A. Strong Near Neighbor Repulsion 

We now introduce disorder in the form of a random site 
dependent chemical potential, /ii — fj,+di where the disor- 
der. Si is uniformly distributed between ±A. A is a tun- 
able parameter characterizing the strength of disorder. 
The first question we want to address is the strength of 
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the disorder necessary to destroy the checkerboard sohd 
phase and what new phase is produced. One can try to 
answer this question with a simple argument based on 
energy balance (Imry-Ma). We start at half filling with 
a perfect checkerboard solid and introduce the site dis- 
order. Suppose that at an empty site there is, due to 
fii, a deep potential well which pulls in a neighboring 
boson. This boson will now have near neighbors which 
it will try push away to rearrange its neighborhood in 
a local checkerboard solid which will, consequently, have 
a mismatch at its boundary with the original checker- 
board. The likelihood of this happening depends on the 
disorder and dimensionality. The energy cost, in d di- 
mensions, due to the mismatch at the boundary scales 
like L"*"^ for a region of length L. On the other hand, 
the energy gained by the bosons by falling into locally fa- 
vorable energy wells scales like L''/^. For d = 1, disorder 
is relevant, no matter how weak it is, it always destroys 
the solid order. For d = 3 or more, the energy cost out- 
weighs the gain and the system maintains checkerboard 
order. The d = 2 case is marginal since both, cost and 
gain, scale like L. Typically, in such marginal cases, the 
conclusion is that disorder will indeed destroy long range 
order but just barely. The correlation length, ^, is very 
long and the system size should be even larger to see the 
effect. 

Numerically, for strong disorder {A/V = 2) we can 
easily see that indeed the gapped checkerboard solid is 
destroyed on lattices as small as L = 12 and is replaced 
by a compressible, glassy insulator with no energy gap. 

Since for very weak disorder, the system size needed 
to see the destruction of solid order is too large for us to 
simulate, we resort to finite size scaling for the interme- 
diate disorder case. Although this does not demonstrate 
directly the validity of the Imry-Ma argument (for which 
very weak disorder is needed) we believe that the results 
we will present are qualitatively similar to what happens 
in the very weak disorder case. 

In Figure ^ we show the density, p, as a function of the 
chemical potential, = (/ii), for i = 8, 12, 14 (/? = 14) 
and L = 20 {(3 = 20) and A/V = 1. For L = 12, 14 and 
20, we average 100 disorder realizations, for i = 8 we 
did 400. We see that the incompressible gapped region 
(k = dp/diJ, = 0) gets smaller as L increases but does not 
quite reach zero. The inset shows the average gap size 
versus L~^. The average gap size was obtained by cal- 
culating the p, fj, curve for each realization, which yields 
the gap size per realization which we then average. An 
alternative method (which washes out important statis- 
tical information) is to calculate the average p, fi curve 
and use it to calculate the average gap. We see clearly 
that the gap tends to zero for a finite, but large, L. This 
suggests that, for these values of V and A, the gap will 
disappear by L « 30. Since L should be greater than 
^ to observe the destruction of the checkerboard order, 
we estimate from this that ^ ~ 30. In Figure we show 
S{7T, tt) as a function of L~^. This, again, shows that the 
checkerboard order is destroyed for systems larger than 
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FIG. 2: p versus p for A/V — 1 and different systems sizes 
showing the shrinking gap. Inset: The average gap versus 
L'\ /3 = 14 for L = 8, 12, 14 and /3 = 20 for L = 20 
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FIG. 3: S'(7r,7r) versus 1/L for A/V = 1. /3 = 14 for L = 
8, 10, 12, 14 and /3 = 20 for L — 20. Long range order seems 
to disappear for L > 30. 



about L = 30. 

To elaborate this further, we show in Fig. ^ the distri- 
butions of the gap sizes for different disorder realizations 
for L ~ 8, 20. We see that for L — 8, the distribution 
is quite narrow and peaked at a nonzero value. How- 
ever, for L = 20, the distribution is very wide and in fact 
peaked at zero indicating that the most probable value 
for the gap is zero. In this case, it is incomplete to discuss 
the "average" of the gap, which is still non-zero. 

In order to characterize further the compressible phase 
which replaces the checkerboard solid, we study the be- 
havior of the Green's function, both for equal and un- 
equal imaginary times. For example, in Figure^ we show 
the equal time Green's function for V — 4.5, L = 10 and 
p ~ 0.56. We see that the Green's function goes to zero 
and is very well fit by an exponential (in fact a hyperbolic 
cosine to account for the periodic boundary conditions). 
This is further evidence that there is no superfluid in this 
phase, it is an insulator. 
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gap gap 

FIG. 4: Gap distribution for different realizations for L — 8 
{f3 = 14, 400 realizations) (a) and L = 20 (/3 = 20, 100 
realizations) (b). V = A = 4.5. As lattice size increases, the 
distribution gets very wide and peaked at 0. 
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FIG. 5: Equal time Green's function as a function of distance 
for L = 10, /3 = 20, V = 4.5, and A/V = 1. The solid line is 
a fit of the form: y = Ao(exp(-(L/2 - x)/Ai) + exp((L/2 - 
x)/Ai)) with Ao = 6.4 x lO""'"' and Ai = 0.552. 



The glassy nature of this phase can be seen in the time- 
dependent Green's function, Ejq. ||. In the Bose glass 
phase, this quantity is predictedElto decay as G(r) ^ 1/r, 
which has been verified numerically in a very different 
contexts. Figure || shows that is also true in this case. 
Therefore, the new phase replacing the gapped checker- 
board solid is an ungapped insulating Bose glass phase. 

Clearly, it is very difficult to examine these issues 
with smaller couplings and disorder: The correlation 
length will be even longer and much larger sizes would 
be needed. However, for moderate disorder, the above 
results demonstrate that, whereas it might appear on a 
finite lattice that the gapped solid phase is still present, 
finite size scaling clearly shows the gap to disappear on 
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FIG. 6: The Green's function, G(r) as a function of imaginary 
time separation, same parameters as Fig. 5. The solid line 
is a fit of the form: G{t) = Aq^t^^ + {[3/2 - t)^^) + A2 
with A(i = 0.04, Ai = -1.0 and A2 = -0.0065. With a 
two parameter fit (excluding Az) we also get a good fit with 
Ao = 0.03 andAi = -1.16. 

large enough systems. We may conclude from this that 
the Imry-Ma argument holds and that disorder, no mat- 
ter how weak, will produce a glassy, compressible, un- 
gapped insulating phase at strong near neighbor cou- 
plings. 

B. Weak Near Neighbor Repulsion 

We now consider the question of what happens when 
the near neighbor repulsion is decreased and only the 
hard core and disorder interactions remain. For soft core 
bosons with no nn interaction, a re-entrant behavior was 
obser'ued for as, a function of the contact repulsion both 
in onea and twocl dimensions. In other words, for fixed 
disorder strength, as the onsite repulsion is increased 
from zero, the superfluid density is at first zero, then 
at some intermediate value of Vo/t the bosons delocal- 
ize and ps takes on a finite value, then for large enough 
V^/t {i.e. approaching the hard core limit) the bosons 
are localized again. In one dimension this happens for 
any amcmnt of disorder, but in two dimensions was only 
reporteda at A/i = 6. This was taken as confirmation 
of the phase diagram presented in reference ^ where the 
bosons were argued to be always localized, even by weak 
disorder, when Vq/I is very large. p. 

The phase diagram in the (t/Vb, /i/Vb) planca should 
however be interpreted with care especially if we want to 
consider the hard core limit. It is a phase diagram at 
constant finite disorder A/Vq and thus I/Vq — > means 
that A/t — > 00, and any arbitrarily weak disorder in units 
of Vb becomes infinitely strong in the hard core limit. 

In fact, figure shows that hard core bosons behave dif- 
ferently, when a finite disorder A/t is considered. While 
the bosons are localized by weak disorder for large nn 
coupling, Vi (and densities that are not very small), as 
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FIG. 7: The Green's function versus r for F = 0, A = 5, /? = 
24, p = 0.1 and L = 8, 10, 12, 16 showing that the system 
is superfluid. The inset shows ps/p versus V for a 12X12 
system. A/t = 5,/3 = 24, p = 0.1 and 250 realizations. Ps/p 
vanishes only for large V . The particle density is p « 0.5. 



this coupling is reduced, the bosons become superfluid 
and stay that way even at Vi = 0. We verified this at 
several particle densities. In other words, the naive ex- 
pectation based on the soft core phase diagram as pre- 
sented in reference ^ at very large onsite repulsion is not 
fulfilled. The soft core phase diagram must be interpreted 
appropriately. 

As we saw for the gap, in the presence of disorder it is 
not enough to consider only average quantities, it is also 
instructive to consider their distribution. In Figure || we 
show the histogram of Psj p for the 250 realizations for 
i = 16, y = 0, A 5i and p « 0.1. We see that for this 
size the distribution is well centered around Psj p ~ 0.18, 
with the same behavior observed for smaller sizes: The 
distribution does not widen nor does the peak tend to 
zero. Therefore, even for this large disorder, the bosons 
are still not localized for F = 0. 

To show that the bosons can be localized at F = if 
the disorder exceeds a threshold value, we show in Fig. ^ 
finite size scaling for Psj p for two cases of disorder, A — 
5t, and 6t. The case of A — 5t was done with p = 
0.1 and clearly shows the system to be superfluid as as 
i — > oo. This is also shown in the equal time Green's 
function. Fig. On the other hand for A = 6t we did 
the simulation at p « 0.73 and we flnd that both Ps/p 
and the Green's function vanish as L — > oo. We also 
verifled this for p « 0.5 and 0.1. We conclude that the 
disorder has to exceed a threshold value (« 6t) in order to 
localize the bosons for weak nn repulsion. Clearly, in the 
no-hopping limit, i — > 0, the critical value of the disorder 
vanishes in agreement with reference |^. 

When the disorder is strong enough to localize the 
bosons at weak coupling and consequently also at strong 
coupling (since even very weak disorder can accomplish 
that), the question arises as to whether there.is re-entrant 
behavior, as observed for soft core bosonsoQ. In other 




FIG. 8; Histogram of ps/p for the 250 reahzations for L 
16, V = 0, A/f = 5, /3 = 24 and p f» 0.1. 
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FIG. 9: Finite size scaling of the superfluid density fraction 
Ps/p for the pure hard-core case = for two different dis- 
order strengths A — 5t and A = 6t, /3 = 20. 
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FIG. 10: Finite size scaling of the equal time Green's function 
at the largest distance for two different disorder strengths 
A = 5t and A = 6t, = 20 
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words, will the bosons go from being localized at weak 
coupling to being delocalized at intermediate values and 
then relocalize at strong coupling? A finite size scaling 
analysis shows that this is not the case: For the hard core 
system, when the disorder is large enough to localize the 
bosons at weak coupling, they will stay localized at all 
couplings. 

This means that the {^/V,t/V) phase diagram is as 
follows. For weak disorder the checkerboard solid is com- 
pletely destroyed and replaced by a compressible insulat- 
ing Bose glass phase. This is easy to understand by ap- 
plying to the two-sublattice structure of the checkerboard 
solid a simple Imry-Ma argument. However, away from 
half filling, there is no such structure and the Imry-Ma 
argument no longer holds. So, with weak disorder, there 
will be a superfluid even for large nn repulsion, when the 
density is far from 0.5. In addition, there is a superfiuid 
for all fillings (except full filling) when the nn repulsion 
is small. For strong disorder, the superfiuid phase ev- 
erywhere is destroyed and replaced by the compressible 
Bose glass phase. 

IV. CONCLUSIONS AND DISCUSSION 

We have used the SSE algorithm to simulate the dis- 
ordered hard core bosonic Hubbard model in two dimen- 
sions to try to understand the interplay and competition 
between interaction and disorder. 

Using finite size scaling, we found that the checker- 
board solid present at strong nn interactions for half fill- 
ing is completely destroyed by any amount of disorder. 
This agrees with the simple Imry-Ma energy balance ar- 
gument. 

Surprisingly, finite size scaling showed that at very 
weak, even vanishing, nn coupling, weak (even intermedi- 



ate) disorder does not localize the hard core bosons. The 
disorder strength, A/t, must be at least of the order of 
6i before the bosons are localized. We found this some- 
what surprising because soft core bosons were argued to 
be always localizedH, even by weak disorder, in the limit 
t/yo ^ oo which is often considered to be equivalent to 
the hard core case. The soft core phase diagranJ3 should 
be interpreted carefully, as discussed in section IIIB. 

Numerically, the bosons were shown to be localized but 
only one value of the disorder was presentedD, A/i = 6, 
which is large. Our numerical results here agree with 
reference^, but show no localization for weaker disorder. 

Another interesting analogy to make is with fermions. 
In one dimension, weak disorder localizes both fermions 
and hard core bosons, which is not surprising in view 
of the fact that they are equivalent in this case. In two 
dimensions, however, this equivalency is lost, and the re- 
sponse to disorder is different. Fermions are (marginally) 
localized by weak disorder while hard core bosons are not. 

Since the hard core bosonic Hubbard model is equiva- 
lent to the spin— i quantum Heisenberg model, the above 
results also hold for this model in the presence of a ran- 
dom external magnetic field. 
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